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Abstract 

The paper develops methods to construct a one-stage optimal design of dilution exper- 
iments under the total available volume constraint typical for bio-medical applications. 
We consider various design criteria based on the Fisher information both is Bayesian 
and non-Bayasian settings and show that the optimal design is typically one-atomic 
meaning that all the dilutions should be of the same size. The main tool is variational 
analysis of functions of a measure which also allows inclusion of other constraints like 
the cost associated with the experiment and its outcome. 
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1 Introduction 

The paper studies a new class of statistical experiments arising, in particular, 
in stem cells research, a very active area of experimental biology. Stem cells 
are the cells produced during early stages of embryonic development and they 
are having capacity to turn into various types of tissue cells. This potentially 
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opens new ways to cure many diseases and explains the huge importance of 



the stem cells research, see, e. g., Mayhall et al. (2004) and the references 



therein. We aim to characterise an optimal design of dilution-type experiments 
under the contraint of the available solution which is typical for experiments 
involving counting stem cells. Specifically, this study originates in studies of 
Hematopoietic or blood stem cells (HSCs) which are the stem cells giving rise 
to all red and white blood cells and platelets. HSCs develop in a mammal's 
embryo in early days from normal cells. There is not much known on the 
biological mechanism which triggers a normal cell to develop into an HSC and 
there is no direct way so far to observe such cells, called pre-HSCs, which soon 
to become HSCs. Pre-HSCs are mostly produced in aorta-gonad-mesonephros 
(ACM) region of the embryonic mesoderm and also in the yolk sac, then colonise 
the liver. A challenging problem is how to detect the pre-HSCs and how many 
of them are present at a different embryo ages. 

In order to estimate the number of pre-HSCs present in a given region, 
experiments on laboratory mice have been conducted which use the following 
signature property of stem cells. A mature HSC is capable to cure a mouse which 
receives a controlled dose of radiation if injected in its blood. A mouse recovers 
(repopulates) , if and only if, it has received at least one HSC in the injected 
dos^J Thus the number of pre-HSCs can be estimated by the so-called limiting 
dilution method: controlled dozes of a substrate containing samples from the 
ACM are injected into irradiated mice and then the number of repopulated 
mice infers on the number of HSCs which developed from the initially present 
pre-HSCs. 

The dilution method has been in the arsenal of biologists for almost a cen- 
tury, at least since McCrady used it for quantitative determination of Bacillus 
coli in water in 1915, |McCrady (1915) Since then many studies used it to esti- 



mate the number of objects of interest in a medium without their direct count, 
Fisher (1922")| |Cochran (1950)} |Ridout (1995)] to name a few. 



So far, several studies have produced estimates for the number of pre-HSCs 
in AGM by using the dilution method which varies between just a few to, per- 



haps, as many as 200, see, e. g., 


Kumaravelu et al. (2002) Gekas et al. (2005) 


Ottersbach and Dzierzak (2005) | 


Bonnefoix and Callanan (2010) and 


Medvinsky 


et al. (2011) Most of the reports tend to focus on the modelling of experimental 



data and on the estimation methods. However, there has been little discussions 
on how to design the experiment to capture the most informative sample. In- 
deed, the experiment would be spoilt, if all the mice repopulate or if all do not. 
The aim of this paper is to find an optimal design of the dilution experiment 
to estimate the mean number of pre-HSCs. To this end, we are going to apply 
recently developed methods of constrained optimisation of functionals of mea- 
sures and the corresponding steepest descent type algorithms for computations. 
It should be stressed that our methodology is generic in the sense that it can be 
applied to other dilution experiments and not only in the stem cells research. 



1 This assumption is rather questionable, so the biologists cautiously speak of one repopu- 
lation unit for this unknown minimal number of HSCs sufficient to cure an irradiated mouse. 
In this study we basically loosely speak of one HSC as of one repopulation unit. 
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Moreover, additional contraints can be incorporated into the model which would 
allow, for instance, to take into account the cost of mice or other materials used 
in the experiment. 

The paper is organised in the following way. Section [2] introduces the dilu- 
tion experiment we are dealing with, description of the corresponding statistical 
model followed by its assumptions, and at the end the optimality criterion func- 
tions. Section [3] lays out the theoretical basis for the optimisation methods we 
are using given that the goal functions are represented as functions of a measure. 
Consecutive sections provide the optimal design of the dilution experiments un- 
der various conditions and various goal functions: in non-Bayesian setting in 
Section |4.1| then under Uniform and Gamma Bayesian priors in Sectio n |4.2| 
We conclude by discussion of our findings and their extensions in Section [5j 

2 Dilution experiment and statistical model 

Description of experiment The dilution experiment on estimation of the 
number of HSCs we address here involved the AGM region of an 11 days old 
mouse embryo. More exactly, in order to make a study more representablc 
and not depending on features of a particular embryo, an engineered AGM 
is made from several such embryos. A substrate of volume V is then made 
from this engineered AGM and the content is thoroughly mixed. Next, n doses 
containing proportions X\, . . . ,x n of the whole V are extracted from all or part 
of the substrate and left for a few days so that pre-HSCs in these doses, if any, 
mature into HSCs. Finally, these n dose^j are injected into n irradiated mice 
(the number n of mice was 30 in this experiment) so that each mouse receives its 
own dose and the mice are put to rest for a few weeks for the doses to take effect. 
After this, the number of repopulated mice, which are the ones having received 
a doze with at least one pre-HSC, is counted and the inference is drawn on the 
total number of pre-HSCs initially present in the AGM. The main question to 
address when designing such an experiment is doses of which volumes should 
be used for available number of mice in order to get best possible quality of 
the statistical estimator? Different criteria could be considered to quantify the 
quality of the estimator. We will consider below most common ones based on 
the Fisher information. 

Statistical model Given the description of the experiment above, the follow- 
ing assumptions can be made. 

• Spatial homogeneity: pre-HSCs were distributed uniformly throughout the 
substrate when the doses were taken. Thus there is no tendency for pairs 
or groups of pre-HSCs either to cluster or to reject one another. This is 
implied by the fact that the substrate was thoroughly mixed just before 
the doses are taken. 

2 These are further diluted to a standard volume, but this, obviously, does not change the 
number of HSCs present before the dilution. 
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• Orderliness: the probability that there are more than one pre-HSCs in a 
small volume dv of substrate has order o(dv). This is a natural assumption 
given that an 11 days old AGM contains about 300 thousand cells and only 
no more than 200 of these are pre-HSCs. 

• Independence: each cell in substrate has the same (small) probability to 
turn into a pre-HSC independently of the other cells. 

• Only a pre-HSC can develop into a mature HSC. So that a dose contains 
at least one HSC at the time of injection to a mouse if and only if there 
was a pre-HSC present in the dose at the time of its extraction from the 
whole substrate. 

• Finally each dose when injected into an irradiated mouse is certain to 
exhibit a positive result (repopulated mouse), whenever the dose contains 
at least one HSC. 

The first three assumptions above suggest that the locations of pre-HSCs in 
the substrate V arc given by a homogeneous Poisson point process. This follows 
from the Poisson limit theorem for thinned point processes, see, e. g., |Daley| 



and Vere- Jones (2008, Sec. 11.3) Indeed, in every subset of volume x of the 
substrate, the number of cells turned into pre-HSCs is close to Poisson distri- 
bution with the parameter proportional to mean number of pre-HSCs in this 
volume, so with Aa; for some A. And because of independence assumption, these 
numbers are independent for disjoint subsets. Changing the units if necessary, 
we assume from now on that the volume of the substrate is 1. The parameter 
A is then the unknown density of the Poisson point process which is also the 
mean number of pre-HSCs in the substrate. Thus we operate with a measur- 
able space carrying point configurations inside a set V C K 3 of volume 1 (the 
space of finite counting measures to onV with the minimal cr-algebra making the 
mappings uj \- > lj(B) measurable for all Borcl B C V) supplied with probability 
distribution Pa so that uj under P\ is a homogeneous Poisson point process 
with density A. 

The doses taken can now be associated with disjoint subsets V\,...,V n of 
V with volumes x,-, i = 1, ...,n. The corresponding numbers of pre-HSCs 
w(Vi), i = 1, . . . , n in the doses are then independent Poisson distributed ran- 
dom variables with parameters Xxi while the total number of pre-HSCs w(V) is 
Poisson distributed with parameter A. 

In the simplest case all the doses have the same volume x < 1/n. The 
probability that a doze is sterile, i.e. it does not contain a pre-HSC, is then 

P = P(u J (V l )=0) = e- x *. (1) 

Thus the total number of non-repopulated mice follows Binomial distribution 
with parameters n,p and the maximum likelihood estimate for the average num- 
ber of HSC A is given by 

A = -^, (2) 
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where p is the proportion of non-repopulated mice provided it is not 0. However, 
the doses need not be necessarily all equal for an optimal design. 

Let Xi {i — 1, n ) be an indicator that a mouse, which received the ith dose 
of volume Xi, has not repopulated. Thus, Xi is & Bernoulli random variable 

Xj |A^Bcrn( e - A ^) (3) 

with the parameter equal to the probability of the ith dose to be sterile. 

Hence, the log-likelihood function for the sequence \ — Oti> ■■■iXn) °f non- 
repopulated and repopulated mice is given by 

n n 

£( X \X, *) = -J2 X* Xx i + E(! " ^) log(l - e"**' ), (4) 

i=l i=l 

where x = (x\, ...,x n ). Maximisation of this expression over A for an observed 
sample x provides a maximum likelihood estimator (MLE) of A for a given 
design x. Our goal here is to determine the optimal design in terms of the doses 
{xi}, according to a suitably chosen optimality criterion, which we describe 
next. 



Optimality criterion functions Recall that for a statistical model which 
depends on a one-dimensional parameter A, the Fisher information is defined 
as 

-d 2 e( x \\,x)' 



J(x; A) 



dX 2 



(5) 



The expectation is taken here w.r.t the random vector x- 
Derived from Q, we have in our case 



/(*;A) = £ r 



(6) 



The Fisher information measures the amount of information that an observable 
sample carries about the unknown parameter, which the likelihood function 
depends upon. On the other hand, it is the inverse of MLE's variance, see, 
e. g., |Everitt (2002)] Thus, maximising the information corresponds to min- 
imising the variance of the MLE. Therefore, maximising the Fisher information 
(|6j) over x under constraint X)"=i x i — •"■ ^ s a useful design criterion. 

It is typical in statistical experiment planning to describe design in term of 
a probabilistic design measure. Typically, the design measure is atomic, so it 
has a form . qjS Xj , where 5 X is the unit mass measure concentrated on a point 
{x}. The design measure reflects the (asymptotic when n — > oo) frequency qi of 
occurrence of the value Xi in the design, see, e. g., Atkinson and Donev (1992) 



By this reason, we will also describe the doses by a measure [i(dx) living on 
[0, 1], albeit not renormalised to have mass 1. Namely, /i = Y]^ mj5 Xj means the 
design when a dose of volume Xj is repeated rrij times. Obviously, we have that 
the total mass constraint JZj—i m j = / n{dx) = n and that we cannot extract 
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more doses than the total volume of the substrate: J2j=i m j x j = / xfi(dx) < 1. 
All integrals here and below are taken over [0, 1], unless specified differently. 
Now the basic optimisation problem for the design of our experiment is 

/e~ Xx 
l _ e -xx x2 ^{dx) -> sup (7) 

over measures \i with supp \x C [0, 1] satisfying 

m([o,i]) = »; (8) 

A design measure describes asymptotic frequencies, so qjn = rrij for a given 
finite n may not be all integers. In this case it is reasonable to consider the 
nearest measure with all rrij £ Z + as an approximation to the optimal solution. 
Or, if necessary, a choice of the optimal measure among such measures can be 
done by evaluation of the goal function at just a few closest approximations of 
this kind to the optimal design measure. 



Bayesian setting Sometimes, there is an additional information available on 
the plausible values of the parameter A which is given in a form of a prior 



distribution Q(dX), see, e. g., Marin and Robert (2007, Sec. 2.2.2) In this case 



the optimality criterion involves taking the expectation of the goal functions 
w.r.t the distribution Q. 

The following criterion functions are typically used for estimation of a single 
parameter, see, e.g., Atkinson and Donev (1992) or Ridout (1995) 



G 2 (m) = E q 7(/x;A) 
G 3 (m)=E q log/(/x; A) 

G 4 (/x) = -E Q (7(/x;A))" 



(10) 
(11) 
(12) 



Criterion function G2 and G3 maximize, under the same constraints (lU and 
([9| , the expectation of the Fisher information and of its logarithm, respectively, 



as used in e.g., in Zacks (1977) and Chaloner and Larntz (1989) The crite 



rion function G4 minimises the expected asymptotic variance of the maximum 
likelihood estimator. 

In the next section will describe a general framework of optimisation of 
functionals of measures and the recently developed techniques for solving such 
optimisation problems. Apart from optimal design of statistical experiments 



Pukelsheim (1983) , Atkinson and Donev (1992) these are frequent in different 
subjects, like spline approximation of curves and geometrical bodies where the 
measure describes the positions of spline points Schneider (1988) maximisation 



the area covered by random geometric objects with the distribution determined 
by a measure Hall (1988)] stochastic search, where a measure determines the 
search strategy Wynn and Zhigljavsky (1994) Zhigljavsky (1991) 



G 



3 Optimisation of functionals of measures 



In this section we summarise necessary information about measures and varia- 
tional analysis on them. Further details on measure theory can be found, e. g. 
in |Dunford and Schwartz (2009)| or |Hille and Philip (1957)] 



Let X be a locally compact separable topological space with the Borel a- 
algebra B of its subsets. Let M (M+) denote the set of signed (respectively, 
non-negative) finite measures on B, i. e. countably additive functions from B 
to K (R+, respectively). M becomes a linear space if the sum of measures 
and multiplication by a number is defined by (77 + v){B) — r/(B) + v{B) and 
{trf){B) = t-q(B) for any t E R and r], v E M. M + is a cone in M since 
fi + v E M + and tfi E M + whenever //, v € M + and t > 0. The support supp /i 
of a positive measure m, is defined as the complement to the union of all open 
sets of zero /^-measure. Measures are orthogonal if their supports are disjoint. 
Any signed measure 77 can be represented as the difference rj + — rf~ of two 
orthogonal non-negative measures i] + ,7]~ E M + (the Jordan decomposition). 
The set M becomes a Banach space if supplied with the total variation norm 
\\ V \\=r 1 +(X)+r,-(X). 

Optimisation of functions defined on a Banach space, which are commonly 
called functionals, relies on the notions of differentiability. In our case a func- 
tional / : M 1— > E is called strongly or Frechet differentiable at v E M if 

f(v + r))- f{y) = Df(v)[rj\ + o(|M|) as | 0, (13) 

where Df(v) is a bounded linear continuous functional on M called the differ- 
ential. 

When a function is strongly differentiable at v then there also exists a weak 
or directional or Gateaux derivative, i. e. 

]imt- 1 [f(v + tr,)-f(v)}=Df(v)[r ] } (14) 

for any 'direction' ?y E M. 

The differential Df{v) often has an integral form: 

D f(y)[v] = J g{ x ;^)v( dx ) 

X 

for some function g{ • ; v) : X 1— > K which is then called the gradient function to 
/ at v. Not all linear functionals are integrals, unless the space X is a finite set 
in which case M can just be identified with an Euclidean space. In most applica- 
tions, however, including the experimental design, differentiable functionals do 
possess a gradient function, so this assumption is not too restrictive in practice. 

In this study we are interested in optimisation of functionals of positive 
measures. A general constrained optimisation problem can be written as follows: 

f(n) — > inf, n E A, (15) 
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where A C M + is a set of measures describing the constraints. If / is strongly 
differentiable then the first order necessary optimality condition states that if 



/i* provides a local minimum in the problem (15) then 



where 



Df(jf)[rj\>0 for all n G T A (/**)> 



Ta(u) = liminf 

yrJ W t 



(16) 



(17) 



is the tangent cone to A at fi. Here the '+' (respectively '— ') operation on 
sets indicates all pairwise sums of (respectively difference between) the points 
from the corresponding sets. The tangent cone is the closure of all admissible 
directions rj £ M at pt meaning that fi+tn £ A for all sufficiently small t > 0, see, 
e. g., |Cominetti (1990)| In other words, derivative in all admissible directions 
should be non-negative at a point of local minimum. For any A of interest, one 
needs to characterise the tangent cone Ta(/x). 

General optimisation theory for functionals of measures has been developed 
in a series of papers Molchanov and Zuyev (2000a) | Molchanov and Zuyev 
(2000b) and Molchanov and Zuyev (2004) | For us here an optimisation with 
finite number of equality and inequality constraints is relevant. 

Consider the following optimisation problem: 



/(//) — » inf, jU G M + subject to 

jH t (m) =0 i = 1, ...,k, k<d 
1-ffj(m)<0 j = k + 1, d. 



i — y 



(18) 
(19) 

l d are strongly diffcrcn- 



where / : M + i-> K and H = (H\, . . . , H d ) : 
tiable functions. Alternatively, the constraints ( 19 ) can be written in the form 
iJ(/i) € C, where C C R d is the cone {y e R^Tyi = • • • = y k = 0, y k+1 < 
0, . • • , Vd < 0}. 



Definition 1 ( |Robinson (1976)[ ). A measure /i is called regular for optimisation 
problem ( 18 ) if the origin of M. d belongs to the interior of the set 



H(v) -C + DH{n)[M + - (i\ C 



(20) 



Robinson's regularity condition which as shown in Zowe and Kurcyusz (1979)] 
guarantees the existence and boundedness of the KuhnTucker vector appearing 
in the next theorem. See also Maurer and Zowe (1979) for the discussion of 



different forms of regularity condition. 

The following theorem gives the first-order necessary conditions for a mini- 



mum in the problem ( 18 ) 



Molchanov and Zuyev (2000a, Th. 4.1) Let /u* be a regular local 
subject to (|19[). Assume that f and H are Frechet 



Theorem 1. 

minimum of f over M + , 

differentiable at ji* , and there exist the corresponding gradient functions g(x, ji 
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and hi(x, fi*), i — 1, d. Then there exists KuhnTucker vector (m, . . . , u<z) € 
E d such that < (resp. Uj — 0) for those j E {k + 1, d} satisfying 
Hj(n*) = (resp. Hj(fi*) <0), such that 



g(x,u* 



for n* -almost all x, 
for all x £ X. 



(21) 



One can show that in the case of finitely many constraints ( 19 ) satisfy- 
ing (21 1, the regularity condition pOl becomes the so-called Mangasarian- 



Fromowitz constraints qualification (see Cominetti (1990, p. 274)), that is the 
linear independence of the gradient functions h\{ •, /z* ),..., hk( •, A**) and the 
existence of a measure 77 € M such that 



hi{x, fj,*) rj(dx) — for all i — 1, k; 
hj(x, fi*) rj(dx) < for all j 6 {k + 1, d} 



verifying 



0. 



(22) 



(23) 



Without the inequality constraints, condition (23 1, trivially holds for 77 being 
the zero measure. 

The design problems we consider naturally fall in the above described general 
framework of optimisation of functionals defined on finite measures. Theorem [T] 
provides the necessary conditions for optimality of a design. Moreover, it allows 
one to easily incorporate into the model other constraints on the optimal design 
measure, if needed. Constraints Q and ^ correspond to linear functionals 
Hi(n) = J n(dx) — n and ^(a 4 ) = Jx fi(dx) — 1 with the corresponding gradient 
functions hi(x; (j,) = 1 and h2(x;fx) = x. These constraints are regular for any 
fi since (p2| and (23) are satisfied, for instance, for a measure 77 = 5q — S\. We 



therefore have the following important corollary of Theorem [T] which we use in 
the next section. Note that we mostly maximise the goal function so that the 



inequalities in (21) change to opposite. 



Theorem 2. Let fi* be a local maximum of a strongly differentiable at fi* 
function f : M + R possessing a gradient function g(x;/j,*), subject to con- 
traints ^ and Then, there exist constants ui and u 2 > where u 2 > if 
J x/i*(dx) — 1 and u-i — if J xfi*(dx) < 1, such that 

1 *\\ —U1+U2X u* -almost everywhere, , . 

g ( x ^ ) 1 ^ , , „ c Y ( 24 ) 

< U\ + u 2 x lor all x € A. 



4 Construction of optimal design 

In this section we apply the necessary condition for extremum of a functional 
of measures to find optimal designs for a range of goal functions and most 
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common prior distributions in the Bayesian settings. First we assume that 
the parameter A, the mean number of HSCs in the substrate, is known from 
previous experiments, and obtain the optimal, in terms of maximisation of the 
Fisher information, design for each A. 

4.1 Optimal design for a fixed average number of HSCs 

Here we are dealing with optimisation problem Q under constraints Q and 
|9|. The goal function is a linear function of ll, so that its differential is the 
function itself with the gradient function 

g 1 (x;\) = Y ^^x 2 , (25) 
not dependent of u,. Note that gi{x) = \~ 2 r(Ax), where 

r(y) = ^— y V 2 - (26) 

The graph of r is shown on Figure [T] It attains its unique maximum at point 
Vinax w 1.59362 and it is strictly concave on [0 

; Umax] • 

The gradient function is just scaled in both directions function r attaining 
its maximum at point y maa; /A, respectively. It follows from Theorem [2] that if 
ll* is a measure at which G\ attains its maximum, then g\(x; A) < u\ + U2X for 
all x € [0, 1] with u 2 > 0. Moreover, g\(x; A) = U\ + u 2 x for x € supp/i*. But 
this is only possible if u\ + u 2 x is a tangent line to g\ at some x* <G [0, y m ax/^} 
and hence the support of optimal fx* consists of only this point x* . Using {[sj) 
and substituting fx = n5 x into Q and ([9| , we come to the optimisation problem 
of one variable: 

r(\x) —> sup over x G [0,1/n], 

so that x* = 1/n for A < y max n and x* — y max /X, otherwise. Thus we have 
proved the following theorem. 

Theorem 3. The optimal design for the problem ^ under constraints ^ and 
(§ consists in n equal doses of volume 1/n for A < d of volume y max /X 

for A > y max ii, where y ma x ~ 1.59362 is the maximum point of the function r 

This indicates that for those A > 1.59362n we need to sample a proportion 
of the substrate 1.59362n/A, and for those A < 1.59362n we have to take all 
the substrate. Therefore, if a good prior point estimate of A is available, a near 
optimal doses of volume x* can be selected, see Figure [T] 

4.2 Optimal design with prior distribution on A 

Typically researches already have an idea on what are the most likely values 
of A. This can be suggested by previous experiments or by analogy with other 
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similar cases and it is given in the form of a prior distribution Q(dX). The 
optimal experimental design now becomes dependent not only on the choice of 
a criterion, but also on the parameters of Q. Formally, the locations of pre-HSCs 
in the substrate now conform to a Cox (or mixed Poisson) point processes with 
parameter measure Q. 

In this section we carry out optimisation of the three goal functions intro- 
duced above in ( 10 1— ( 12 ) for the most common prior distributions: Uniform and 
Gamma. 



Uniform prior distribution. The Uniform prior distribution is appropriate 
when there is no knowledge on the mean number A of HSCs in substrate apart 
from its range. Certainly, A should be greater than 1 since HSCs do develop 
in the AGM. So, it is reasonable to assume that A ~ Unif(l,u), where u > 1 
is a known upper bound. It was already mentioned that the number of HSCs 
hardly exceeds 200, so one can set u = 200, or, keeping in mind the variance of 
the Poisson distribution, u = 170, for instance. 



The goal function ( 10 ) now becomes 

u 1 

1 f f e 



G 2 (p;u) = - — - / / ^z^x n(dx)d\. 



u — 1 J J 1 — e 
1 o 



By Fubini 's theorem we can change the order of integrals above to arrive at 

l 

1 /" 1 g — UX 

G 2 (/i;«) = r / x log ,~ e _ x n{dx) . (27) 

u — 1 ./ 1 — e x 



Thus G2(p;u) is a linear functional of fi with the gradient function 

x 1 — e~ ux 

92(x; u) = log — . (28) 

u — 1 1 — e x 

This function varies very little for all practically interesting values of u > 5 
with the maximum value attained around 0.69 — 0.7, see Figure [2j Reasoning 
the same way as we did in Section |4.1[ we can conclude that for any n > 2 
the optimal design is attained on one atom design measure concentrated on the 
point 1/n (c/. Plot (c) in Figure[I]). So the whole volume of the substrate should 
be divided into n equal doses for this goal function. 



Consider now optimisation criterion (11 1 which now takes the form 



l 



G 3 (m; u) = f logY I i e _ e _ Xx x-,,uL- i ) ( /A. (2<i) 



1 o 
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This goal function is not linear, but still Frechet differentiable with the gradient 
function given by 

u 

1 f e~ Xx 

g 3 (x,fi;u) = / I~ (n;X)- j-x 2 dX, (30) 

u — 1 J 1 — e Ax 

l 

where 

l ^ 

7( M ;A) = J ^ J_ Xx x 2 n{dx) (31) 
o 

is the Fisher information written in terms of a design measure. Since the gradi- 
ent now depends on fi, numeric methods have to be employed to find the optimal 
design for a given value of the upper bound u and the number of mice n. This 
was done by means of an R-library medea which finds an optimal solution to a 
measure optimisation problem under linear constrains of equality type. Since 
medea docs not yet allow for inequality constraints, to deal with ([9]), the proce- 
dure looks for an optimum measure for a given value b of the integral J xfi(dx) 
and then optimises over b < 1. The R-code is freely available from one of the 
authors web-pag^] 

Numeric experiments conducted for various values of u show that the ob- 
tained optimal solution is always one-atom as in the previous cases. Typically, 
the numeric solution gives two atoms at the neighbouring points of the discre- 
tised space [0, 1] for x which indicates a single atom situated in between these 
grid points, see Figure [3] Although we cannot formally prove that the optimum 
design is one- atom, such designs (i. e. equal doses) are certainly of an interest. 

One-atom measures under constraint (fSp have a form /x = n5 x for some 



x € (0, 1/n]. Making use of (26), we come to a one variable optimisation 
problem: maximise 

Ga(n6 x ) = log 7i — 2 Eq log A + Eq log r (Ax) (32) 

subject to x* < 1/n and Q is Unif[l,it]. Only the last term depends on x, so 
the equivalent problem is to maximise 

u 

G 3 (x;u) =Enlogr(Aa;) = / logr(Aa;)dA. 

u — lj 
l 

As u grows, the maximum x max of this function, which can easily be computed 
numerically, approaches 0, see the right plot on Figure |4j So when x max > 1/n, 
the constrained maximum is attained at the point 1/n. Otherwise, for large 
u, x max < 1/n and the solution is to take n doses of volume x max . This is 
exemplified at Figure [4] for n = 30: the optimal dose is given by 

= (33) 



^http: / /www. math. chalmers.se/-sergei 
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where u* 90.71. This indicates that for those u < u* one needs to take all 
the substrate to make n equal doses, and the volume nx max otherwise. 



Finally, consider objective function G4 given by ( 12 ) which corresponds to 
the expected asymptotic variance of the maximum likelihood estimator. In the 
case of a Uniform prior A ~ U n if (1, it), 

u 

G 4 (^; u) = J j-i(A, fi) dX, (34) 



where /(/i; A) is given by (31). Again, G4 is Frechet differentiable with the 
gradient function 



g i {x,n;u) = / / (A,/i)- -^-x dX (35) 

u — 1 J 1 — e 

1 

and a numeric procedure should be used to find the optimal measure for any 
given values of u and n. Similarly to the case of objective function G3 above, 
numeric experiments show that the optimal measure is one atomic, although we 
cannot show this rigorously. In the class of one-atomic measures \x — nS x the 
goal function simplifies to 

G i {n8 x] u) = - eUX - eX ~ x{u - l) (36) 
nx a (u — I) 

Similarly to G3 above, the point of maximum x max of this function approaches 
when u grows, see the right plot on Figure [5j So that when u is such that 
Xmax > 1/n the optimal dose is 1 jn (the whole substrate is used) , otherwise the 
optimal dose is x rnax . For example, for n = 30 mice, the optimal dose is also 



given by (33) but with u* m 64.48 this time. 



Gamma prior distribution. Gamma prior distribution arises naturally as 
a posterior distribution for the Poisson parameter A which has a Uniform prior 
distribution. Since the probability of having k pre-HSCs in the substrate is 

u 

PMV) = k) = E Q PMU) = k\X) = J dX, (37) 

1 

The posterior p.d.f. for A is then 

&CW-*)- t |/^ ft «.^, (38) 

which is close T(k+1, 1) distribution for large u. So once an estimate for the total 
number of pre-HSCs is available, the Gamma prior would be a reasonable for 



13 



A in the subsequent experiment. Moreover, posterior distribution for a Gamma 
prior r(a, 1) given k pre-HSCs will be T(a + k, 2) and so on. So generally, 
A ~ T(a,j3) with rate parameter f3 S N could be a reasonable assumption. 



Under this assumption, Criterion ( 10 ) becomes 

1 oo 



-Xx 



1-e 



-Xx 



x 2 dX fj,(dx), 



o o 



where 



x «-i e -0\_ 



So Gi is a linear functional with the gradient function 



(39) 



(40) 



^-x 2 f a ^(X)d\ = / r(Xx)X- 2 f a ^(X)dX 



B 2 



(a-l)(a-2) 



r(Xx)f a ~ 2 ,p(X)dX 



(a-l)(a-2) 



r(A)/ Q _ 2i/3/x (A)dA. (41) 



Thus the situation here is similar to the case of the Uniform distribution: de- 
pending on whether the maximum of this function is below or above 1/n one 
needs to take equal doses of volumes equal the point of maximum x max or 1/n. 

Consider as an example the first-iteration case when (3 = 1. The gradient 
function and the point of maximum are shown in Figure [6] The optimal design 
is given by 



a < a* 

$ 

a > a 



(42) 



where a* w 49.68. 



Under Gamma prior (40), criterion function (11) takes the following form: 

(43) 



G 3 (/z;a,/3)= / f a ,fs(X) log/(/i; A)dA, 



where 7(/i; A) is given by (31). Note that G3 is not a linear functional of a 
measure, however it is Frechet differentiable with the gradient function 



g 3 (x,(x;a,0) = I 1 (p;X) 



-Xx 



f a . p (X)dX. 



(44) 
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Thus, as before numeric methods for /3 = 1 and a given values of a and n should 
be employed. Our experiments show that the optimum measure given by the 



steepest descent algorithm still contains only one atom. So, optimising (43 1 for 



one atomic measures nS x over x € (0, 1/n] leads to the design given by (42) 
with a* — 47.73. Figure [7] shows G^{n8 x ]a) and the optimal equal doses for 
different values of a. 

Finally, consider G4 in (fl2|) with Gamma prior distribution on A: 



GtfaaJ) = -j I- 1 (fi;X)f a ^(X)d\, (45) 
o 

which is Frechet differentiable with the gradient function 

oo 

/g-Ax 
j ~ 2 (m; A) 1 - e _ Xx x f at p(\)d\. (46) 



Here also, our numeric experiments for /? = 1 and any given values of a and n, 
show that the optimal measure contains a single atom. Thus, in the class of one 
atomic measures, (45) simplifies to 



Gi(n6 x ;a) 



1 f 1 - e~ Xx 





(47) 



which leads to the one- atom optimal design (42) with a* ps 45.75, see Figure [81 



5 Discussion 

We have considered dilution experiments with volume constraints typical in 
biological and medical research. An important particularity of the experiments 
we consider here is that the time delay necessary for pre-HSCs to develop and 
then for the injected doses to take effect in mice prevents from planing multiple 
stage experiments: the estimation should necessarily be done from the first and 
only stage. As we have seen, in all the considered cases of the goal function, 
the optimal design is attained on a one-point measure, meaning that all the 
doses should have equal volume. This parallels the well known result about the 
D-optimal design measure for a linear regression model: The Kiefer-Wolfowitz 
theorem assures that such measure is atomic with the number of atoms to be 
at most the number of contraints plus one, see, e. g., |Kiefer and W olfowitz 



(1960) In our case we have one volume constraint (in addition to the measure 
to have a fixed total mass) so the number of atoms is at most two. Indeed, the 
gradient functions we observed are convex on the interval from to the point of 
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maximum, so the only possibility to satisfy the necessary optimality condition 
given in Theorem |2j is for the optimal measure to have only one atom. 

In practical terms, if the number of HSC expected to be not very large (a few 
dozens or less), the whole substrate should be used to derive the doses, otherwise 
only part of it. We have characterised above what is 'not very large' and how 
much the doses should be diluted. We have fixed some of the parameters here 
(the number of mice n = 30, (3 = 1), driven by practical applications to stem 
cell research. For other applications and other values of these parameters one 
make use the computer codes the authors make freely available for download. 

Advantage of measure optimisation approach is that additional requirements 
can be easily incorporated into the goal function or added as further constraints, 
e. g. limited cost associated with non-repopulated mice and/or the cost of the 
whole experiment if all mice are repopulated or all are not repopulated. 
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(a) 



(b) 




Figure 1: Plot (a): function r, the maximum is attained at the point y ma x ~ 
1.59362. (b): the optimal dose volume for n = 30 mice as a function of A: 
1/30 for A < 47.8 and 0.05312/A otherwise. The line u% + u-ix and the gradient 
function gi(x;X) satisfying conditions (24) for A = 20 (Plot (c), the tangent 
line at point 1/30) and for A = 100 (Plot (d), the tangent line at the point of 
maximum y ma x/\ < 1/30). 
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Figure 2: Function (u - l)g 2 (x;u) for the values u = 10,20,200,1000 (in the 
order of increase). The last two curves are almost indistinguishable. 
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Figure 3: Numeric solution obtained by medea for u — 120 (left pair of plots) 
and for u — 20 (right pair) for n = 30 mice. For u — 120, one atomic optimal 
measure is 30^., where x* = 0.02522. The numeric solution is two atom 
measure 29.97 So. 025 + 0.03 <5o. 026- I n a similar way, we have for u — 20, one 
atomic optimal measure 30 5 X * with x* = 1/30. Moreover, the numeric solution 
is two atom measure 20.07 So. 033 + 9.93 So. 034- So, the total mass of 30 was 
distributed between the grid points surrounding the true atom position. 
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Figure 4: Left plot: function G3 for u = 50, 100 and 200. As u grows, the point 
of maximum x max approaches 0. The vertical dashed line is through the point 
1/n with n = 30. Right plot: x max as a function of u (horizontal dashed line is 
through 1/n for n = 30). 
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Figure 5: Left plot: function — \og(\nGi(n8 x ; u)\) for u = 20,50 and 100 (from 
upper to lower curve). As u grows, the point of maximum x max approaches 0. 
The vertical dashed line is through the point 1/n with n = 30. Right plot: x max 
as a function of u (horizontal dashed line is through 1/n for n = 30). 



22 




200 



Figure 6: Left plot: function g2(x) for /? = 1 and a = 30, 50, 75, 100 (from top 
to bottom). The vertical dashed line is through the point 1/n with n — 30. 
Right plot: x m ax as a function of a (horizontal dashed line is through 1/n for 
n = 30 intersection the curve at point u* w 49.68). 
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Figure 7: Left plot: function G 3 (nS x ; a) for /3 = 1 and a — 50,100 and 200. 
Similarly to the Uniform prior case as a grows, the point of maximum x max 
approaches 0. The vertical dashed line is through the point 1/n with n = 30. 
Right plot: x max as a function of a (horizontal dashed line is through 1/n for 
n = 30). 



23 




Figure 8: Left plot: function — \og(\nG 4(716 x ; u)\) for /? = 1 and a = 20,50 and 
100. Similar to the Uniform case for G4, as a grows, the point of maximum 
x-max approaches 0. The vertical dashed line is through the point 1/n with 
n = 30. Right plot: x max as a function of a (horizontal dashed line is through 
1/n for n = 30). 
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